home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / pcutils / os2 / rt / quad.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-31  |  11.2 KB  |  450 lines
  1. /*
  2.  
  3. QUAD.C  Arbitrary quadratic logic
  4.  
  5.                2   2   2
  6. Defined solid for :- ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j<=0
  7.  
  8. Ellisoids, spheres, infinitely long cylinders, double cones and planes are
  9. all special case of this equation.
  10.  
  11. */
  12.  
  13. /*...sincludes:0:*/
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <stddef.h>
  17. #include <malloc.h>
  18. #include <memory.h>
  19. #include <math.h>
  20. #include "standard.h"
  21. #include "rt.h"
  22. #include "vector.h"
  23. #include "sil.h"
  24. #define    _QUAD_
  25. #include "quad.h"
  26.  
  27. /*...vrt\46\h:0:*/
  28. /*...vvector\46\h:0:*/
  29. /*...vsil\46\h:0:*/
  30. /*...vquad\46\h:0:*/
  31. /*...e*/
  32.  
  33. /*...screate_quad        \45\ create any general quadratic:0:*/
  34. QUAD    *create_quad(
  35.     double a, double b, double c,
  36.     double d, double e, double f,
  37.     double g, double h, double i,
  38.     double j
  39.     )
  40.     {
  41.     QUAD    *quad;
  42.  
  43.     if ( (quad = malloc(sizeof(QUAD))) == NULL )
  44.         return NULL;
  45.  
  46.     quad->a = a;
  47.     quad->b = b;
  48.     quad->c = c;
  49.     quad->d = d;
  50.     quad->e = e;
  51.     quad->f = f;
  52.     quad->g = g;
  53.     quad->h = h;
  54.     quad->i = i;
  55.     quad->j = j;
  56.  
  57.     return quad;
  58.     }
  59. /*...e*/
  60. /*...screate_ellipsoid   \45\ create an ellipsoid of a given size at the origin:0:*/
  61. QUAD    *create_ellipsoid(double rx, double ry, double rz)
  62.     {
  63.     return create_quad(1.0 / (rx * rx), 1.0 / (ry * ry), 1.0 / (rz * rz),
  64.                0.0, 0.0, 0.0,
  65.                0.0, 0.0, 0.0,
  66.                -1.0);
  67.     }
  68. /*...e*/
  69. /*...screate_x_ell_cyl   \45\ create elliptical cylinder along x axis of given radii:0:*/
  70. QUAD    *create_x_ell_cyl(double ry, double rz)
  71.     {
  72.     return create_quad(0.0, 1.0 / (ry * ry), 1.0 / (rz * rz),
  73.                0.0, 0.0, 0.0,
  74.                0.0, 0.0, 0.0,
  75.                -1.0);
  76.     }
  77. /*...e*/
  78. /*...screate_y_ell_cyl   \45\ create elliptical cylinder along y axis of given radii:0:*/
  79. QUAD    *create_y_ell_cyl(double rx, double rz)
  80.     {
  81.     return create_quad(1.0 / (rx * rx), 0.0, 1.0 / (rz * rz),
  82.                0.0, 0.0, 0.0,
  83.                0.0, 0.0, 0.0,
  84.                -1.0);
  85.     }
  86. /*...e*/
  87. /*...screate_z_ell_cyl   \45\ create elliptical cylinder along z axis of given radii:0:*/
  88. QUAD    *create_z_ell_cyl(double rx, double ry)
  89.     {
  90.     return create_quad(1.0 / (rx * rx), 1.0 / (ry * ry), 0.0,
  91.                0.0, 0.0, 0.0,
  92.                0.0, 0.0, 0.0,
  93.                -1.0);
  94.     }
  95. /*...e*/
  96. /*...screate_x_cyl       \45\ create cylinder along x axis of given radii:0:*/
  97. QUAD    *create_x_cyl(double r)
  98.     {
  99.     return create_x_ell_cyl(r, r);
  100.     }
  101. /*...e*/
  102. /*...screate_y_cyl       \45\ create cylinder along y axis of given radii:0:*/
  103. QUAD    *create_y_cyl(double r)
  104.     {
  105.     return create_y_ell_cyl(r, r);
  106.     }
  107. /*...e*/
  108. /*...screate_z_cyl       \45\ create cylinder along z axis of given radii:0:*/
  109. QUAD    *create_z_cyl(double r)
  110.     {
  111.     return create_z_ell_cyl(r, r);
  112.     }
  113. /*...e*/
  114. /*...screate_x_ell_cone  \45\ create elliptical cone along x axis of given radii:0:*/
  115. QUAD    *create_x_ell_cone(double ky, double kz)
  116.     {
  117.     return create_quad(-1.0, ky * ky, kz * kz,
  118.                0.0, 0.0, 0.0,
  119.                0.0, 0.0, 0.0,
  120.                0.0);
  121.     }
  122. /*...e*/
  123. /*...screate_y_ell_cone  \45\ create elliptical cone along y axis of given radii:0:*/
  124. QUAD    *create_y_ell_cone(double kx, double kz)
  125.     {
  126.     return create_quad(kx * kx, -1.0, kz * kz,
  127.                0.0, 0.0, 0.0,
  128.                0.0, 0.0, 0.0,
  129.                0.0);
  130.     }
  131. /*...e*/
  132. /*...screate_z_ell_cone  \45\ create elliptical cone along z axis of given radii:0:*/
  133. QUAD    *create_z_ell_cone(double kx, double ky)
  134.     {
  135.     return create_quad(kx * kx, ky * ky, -1.0,
  136.                0.0, 0.0, 0.0,
  137.                0.0, 0.0, 0.0,
  138.                0.0);
  139.     }
  140. /*...e*/
  141. /*...screate_x_cone      \45\ create cone along x axis of given radii:0:*/
  142. QUAD    *create_x_cone(double k)
  143.     {
  144.     return create_x_ell_cone(k, k);
  145.     }
  146. /*...e*/
  147. /*...screate_y_cone      \45\ create cone along y axis of given radii:0:*/
  148. QUAD    *create_y_cone(double k)
  149.     {
  150.     return create_y_ell_cone(k, k);
  151.     }
  152. /*...e*/
  153. /*...screate_z_cone      \45\ create cone along z axis of given radii:0:*/
  154. QUAD    *create_z_cone(double k)
  155.     {
  156.     return create_z_ell_cone(k, k);
  157.     }
  158. /*...e*/
  159. /*...scopy_quad          \45\ make a copy of a quad:0:*/
  160. QUAD    *copy_quad(QUAD *quad)
  161.     {
  162.     QUAD    *copy;
  163.  
  164.     if ( (copy = malloc(sizeof(QUAD))) == NULL )
  165.         return NULL;
  166.  
  167.     memcpy(copy, quad, sizeof(QUAD));
  168.     return copy;
  169.     }
  170. /*...e*/
  171. /*...sdestroy_quad       \45\ destroy a quad:0:*/
  172. void    destroy_quad(QUAD *quad)
  173.     {
  174.     free(quad);
  175.     }
  176. /*...e*/
  177.  
  178. /*...strans_x_quad       \45\ translate by amount in x direction:0:*/
  179. void    trans_x_quad(QUAD *quad, double t)
  180.     {
  181.     quad->j += (quad->a * t * t - quad->g * t);
  182.     quad->i -= quad->f * t;
  183.     quad->h -= quad->d * t;
  184.     quad->g -= 2.0 * quad->a * t;
  185.     }
  186. /*...e*/
  187. /*...strans_y_quad       \45\ translate by amount in y direction:0:*/
  188. void    trans_y_quad(QUAD *quad, double t)
  189.     {
  190.     quad->j += (quad->b * t * t - quad->h * t);
  191.     quad->i -= quad->e * t;
  192.     quad->h -= 2.0 * quad->b * t;
  193.     quad->g -= quad->d * t;
  194.     }
  195. /*...e*/
  196. /*...strans_z_quad       \45\ translate by amount in z direction:0:*/
  197. void    trans_z_quad(QUAD *quad, double t)
  198.     {
  199.     quad->j += (quad->c * t * t - quad->i * t);
  200.     quad->i -= 2.0 * quad->c * t;
  201.     quad->h -= quad->e * t;
  202.     quad->g -= quad->f * t;
  203.     }
  204. /*...e*/
  205. /*...sscale_x_quad       \45\ scale by factor in x direction:0:*/
  206. void    scale_x_quad(QUAD *quad, double factor)
  207.     {
  208.     quad->a /= (factor * factor);
  209.     quad->d /= factor;
  210.     quad->f /= factor;
  211.     quad->g /= factor;
  212.     }
  213. /*...e*/
  214. /*...sscale_y_quad       \45\ scale by factor in y direction:0:*/
  215. void    scale_y_quad(QUAD *quad, double factor)
  216.     {
  217.     quad->b /= (factor * factor);
  218.     quad->d /= factor;
  219.     quad->e /= factor;
  220.     quad->h /= factor;
  221.     }
  222. /*...e*/
  223. /*...sscale_z_quad       \45\ scale by factor in z direction:0:*/
  224. void    scale_z_quad(QUAD *quad, double factor)
  225.     {
  226.     quad->c /= (factor * factor);
  227.     quad->e /= factor;
  228.     quad->f /= factor;
  229.     quad->i /= factor;
  230.     }
  231. /*...e*/
  232. /*...srot_x_quad         \45\ rotate about x axis by given angle:0:*/
  233. void    rot_x_quad(QUAD *quad, double angle)
  234.     {
  235.     double    b    = quad->b;
  236.     double    c    = quad->c;
  237.     double    d    = quad->d;
  238.     double    e    = quad->e;
  239.     double    f    = quad->f;
  240.     double    h    = quad->h;
  241.     double    i    = quad->i;
  242.     double    ca   = cos(angle);
  243.     double    sa   = sin(angle);
  244.     double    caca = ca * ca;
  245.     double    sasa = sa * sa;
  246.     double    saca = sa * ca;
  247.  
  248.     quad->b = b * caca + c * sasa - e * saca;
  249.     quad->c = b * sasa + c * caca + e * saca;
  250.     quad->d = d * ca - f * sa;
  251.     quad->e = 2.0 * (b - c) * saca + e * (caca - sasa);
  252.     quad->f = d * sa + f * ca;
  253.     quad->h = h * ca - i * sa;
  254.     quad->i = h * sa + i * ca;
  255.     }
  256. /*...e*/
  257. /*...srot_y_quad         \45\ rotate about y axis by given angle:0:*/
  258. void    rot_y_quad(QUAD *quad, double angle)
  259.     {
  260.     double    a    = quad->a;
  261.     double    c    = quad->c;
  262.     double    d    = quad->d;
  263.     double    e    = quad->e;
  264.     double    f    = quad->f;
  265.     double    g    = quad->g;
  266.     double    i    = quad->i;
  267.     double    ca   = cos(angle);
  268.     double    sa   = sin(angle);
  269.     double    caca = ca * ca;
  270.     double    sasa = sa * sa;
  271.     double    saca = sa * ca;
  272.  
  273.     quad->a = a * caca + c * sasa + f * saca;
  274.     quad->c = a * sasa + c * caca - f * saca;
  275.     quad->d = e * sa + d * ca;
  276.     quad->e = e * ca - d * sa;
  277.     quad->f = 2.0 * (c - a) * saca + f * (caca - sasa);
  278.     quad->g = i * sa + g * ca;
  279.     quad->i = i * ca - g * sa;
  280.     }
  281. /*...e*/
  282. /*...srot_z_quad         \45\ rotate about z axis by given angle:0:*/
  283. void    rot_z_quad(QUAD *quad, double angle)
  284.     {
  285.     double    a    = quad->a;
  286.     double    b    = quad->b;
  287.     double    d    = quad->d;
  288.     double    e    = quad->e;
  289.     double    f    = quad->f;
  290.     double    g    = quad->g;
  291.     double    h    = quad->h;
  292.     double    ca   = cos(angle);
  293.     double    sa   = sin(angle);
  294.     double    caca = ca * ca;
  295.     double    sasa = sa * sa;
  296.     double    saca = sa * ca;
  297.  
  298.     quad->a = a * caca + b * sasa - d * saca;
  299.     quad->b = a * sasa + b * caca + d * saca;
  300.     quad->d = 2.0 * (a - b) * saca + d * (caca - sasa);
  301.     quad->e = f * sa + e * ca;
  302.     quad->f = f * ca - e * sa;
  303.     quad->g = g * ca - h * sa;
  304.     quad->h = g * sa + h * ca;
  305.     }
  306. /*...e*/
  307.  
  308. /*...sisects_reqd_quad   \45\ max number of isects we will generate:0:*/
  309. int    isects_reqd_quad(QUAD *quad)
  310.     {
  311.     quad=quad; /* Suppress 'unref arg' compiler warning */
  312.  
  313.     return 4;
  314.     }
  315. /*...e*/
  316. /*...sintersect_quad     \45\ intersect quadratic:0:*/
  317. /*
  318.  
  319. Those who don't like maths should page down a few times very rapidly!
  320.  
  321. Any point along the line we are interested in is of the form p + tq.
  322.                                  ~    ~
  323.       2   2   2
  324. Given:    ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j=0
  325.  
  326. Subs:    x = x +tx        y = y +ty        z = z +tz
  327.          p   q             p   q             p   q
  328.  
  329.        2   2   2                    2
  330. Gives:    (ax +by +cz +dx y +ey z +fz x )t  +
  331.        q   q   q   q q   q q   q q
  332.  
  333.     (2[ax x +by y +cz z ]+d[x y +x y ]+e[y z +y z ]+f[z x +z x ]+gx +hy +iz )t +
  334.          p q   p q   p q     p q  q p     p q  q p     p q  q p    q   q   q
  335.  
  336.        2   2   2
  337.     (ax +by +cz +dx y +ey z +fz x +gx +hy +iz +j) = 0
  338.        p   p   p   p p   p p   p p   p   p   p
  339.  
  340.  
  341.       2   2   2                   
  342. Opt:    ax +by +cz +dx y +ey z +fz x            costs 12x's and 5+'s
  343.       q   q   q   q q   q q   q q
  344.  
  345.     (ax +dy )x +(by +ez )y +(cz +fx )z        costs 9x's and 5+'s
  346.        q   q  q    q   q  q    q   q  q
  347.  
  348.  
  349.       2   2   2
  350. Opt:    ax +by +cz +dx y +ey z +fz x +gx +hy +iz +j    costs 15x's and 9+'s
  351.       p   p   p   p p   p p   p p   p   p   p
  352.  
  353.     (ax +dy +g)x +(by +ez +h)y +(cz +fx +i)z +j    costs 9x's and 9+'s
  354.        p   p    p    p   p    p    p   p    p
  355. */
  356.  
  357. /*...svalue_of_quad:0:*/
  358. /*
  359. Use this to determine value of quadratic at given point in space.
  360.  
  361.           2   2   2
  362. Opt:    ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j        costs 15x's and 9+'s
  363.  
  364.     (ax+dy+g)x+(by+ez+h)y+(cz+fx+i)z+j        costs 9x's and 9+'s
  365.  
  366. */
  367.  
  368. static double value_of_quad(void *shapeinfo, VECTOR v)
  369.     {
  370.     QUAD    *quad = (QUAD *) shapeinfo;
  371.     double    x = v.x, y = v.y, z = v.z;
  372.  
  373.     return (quad->a * x + quad->d * y + quad->g) * x +
  374.            (quad->b * y + quad->e * z + quad->h) * y +
  375.            (quad->c * z + quad->f * x + quad->i) * z +
  376.             quad->j;
  377.     }
  378. /*...e*/
  379.  
  380. void    intersect_quad(QUAD *quad, VECTOR p, VECTOR q, SIL *sil)
  381.     {
  382.     double    a  = quad->a;
  383.     double    b  = quad->b;
  384.     double    c  = quad->c;
  385.     double    d  = quad->d;
  386.     double    e  = quad->e;
  387.     double    f  = quad->f;
  388.     double    g  = quad->g;
  389.     double    h  = quad->h;
  390.     double    i  = quad->i;
  391.     double    j  = quad->j;
  392.     double    qa = (a * q.x + d * q.y) * q.x +
  393.              (b * q.y + e * q.z) * q.y +
  394.              (c * q.z + f * q.x) * q.z;
  395.     double    xx = (a * p.x * q.x + b * p.y * q.y + c * p.z * q.z);
  396.     double    qb = xx + xx +
  397.              d * (p.x * q.y + q.x * p.y) +
  398.              e * (p.y * q.z + q.y * p.z) +
  399.              f * (p.z * q.x + q.z * p.x) +
  400.              g * q.x + h * q.y + i * q.z;
  401.     double    qc = (a * p.x + d * p.y + g) * p.x +
  402.              (b * p.y + e * p.z + h) * p.y +
  403.              (c * p.z + f * p.x + i) * p.z + j;
  404.  
  405.     intersect_quadratic_t(qa, qb, qc, p, q,
  406.                   (void *) quad, value_of_quad, sil);
  407.     }
  408. /*...e*/
  409. /*...snormal_to_quad     \45\ find normal of surface at a given point:0:*/
  410. /*
  411.  
  412. Use partial derivatives to find normal at a given point.
  413.  
  414. Given:      2   2   2
  415.     ax +by +cz +dxy+eyz+fzx+gx+hy+iz+j=0
  416.  
  417. d/dx:    2ax+dy+fz+g
  418.  
  419. d/dy:    2by+dx+ez+h
  420.  
  421. d/dz:    2cz+ey+fx+i
  422.  
  423. */
  424.  
  425. VECTOR    normal_to_quad(QUAD *quad, VECTOR p)
  426.     {
  427.     double    a = quad->a;
  428.     double    b = quad->b;
  429.     double    c = quad->c;
  430.     double    d = quad->d;
  431.     double    e = quad->e;
  432.     double    f = quad->f;
  433.     double    g = quad->g;
  434.     double    h = quad->h;
  435.     double    i = quad->i;
  436.     double    x = p.x;
  437.     double    y = p.y;
  438.     double    z = p.z;
  439.     double    ax = a * x;
  440.     double    by = b * y;
  441.     double    cz = c * z;
  442.     VECTOR    normal;
  443.  
  444.     normal.x = ax + ax + d * y + f * z + g;
  445.     normal.y = by + by + d * x + e * z + h;
  446.     normal.z = cz + cz + e * y + f * x + i;
  447.     return normal;
  448.     }
  449. /*...e*/
  450.